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Abstract: In this paper, we show how to recover the motion of an observer relative 
to a planar surface directly from image brightness derivatives. We do not compute the 
optical flow as an intermediate step. We derive a set of nine non-linear equations using 
a least-squares formulation. A simple iterative scheme allows us to find either of two 
possible solutions of these equations. An initial pass over the relevant image region is 
used to accumulate a number of moments of the image brightness derivatives. All of the 
quantities used in the iteration can be efficiently computed from these totals, without the 
need to refer back to the image. A new, compact notation allows us to show easily that 
there are at most two planar solutions. 
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1. Introduction 

The problem of determining rigid body motion and surface structure from image data has 
been the topic of many recent research papers in the area of machine vision. Two types 
of approaches, discrete and continuous, have been pursued. In the discrete approach, 
information about a finite number of points is used to reconstruct the motion [4,5,9,10,13- 
lb]. To do this, one has to identify and match feature points in a sequence of images. The 
minimum number of points required depends on the number of images. In the continuous 
approach, the optical flow that is the apparent velocity of brightness patterns is used at 
every image point [1-3,12,17]. In general, the computation of the local flow field involves 
not only a least-squares formulation that uses a constraint equation between the optical 
flow and the image brightness derivatives at every image point, but also an additional 
constraint, such as the smoothness of the flow field [7,8]. The additional constraint is 
required because information is lost when the 3-D world is projected onto a 2-D image. 
Locally, the brightness variations in time-varying images only provide one constraint 
on the components of the optical flow. While the continuous approach requires more 
computation, it is robust with respect to errors in the optical flow data. 

Much of the research work in recovering surface structure and motion has been re¬ 
stricted to using the optical flow information in the image. In this paper, we determine 
the motion of an observer relative to a planar surface directly from the image brightness 
derivatives, without having to compute the optical flow as an intermediate step. We 
restrict ourselves to planar surfaces since only three parameters are needed to specify 
the surface structure. We will be considering less restricting constraints on the surface 
structure in a later paper. 

2. Preliminaries 

We first recall some details about perspective projection, the motion field, the brightness 
change constraint equation, rigid body motion and planar surfaces. This we do using 
vector notation in order to keep the resulting equations as compact as possible. 

2.1. Perspective Projection 

Let the center of projection be at the origin of a Cartesian coordinate system. Without 
loss of generality we assume that the effective focal length is unity. The image is formed 
on the plane z — 1, parallel to the xy-plane, that is, the optical axis lies along the z-axis. 
Let R be a point in the scene. Its projection in the image is r, where 



The z-component of r is clearly equal to one, that is r • z = 1. 


2.2. Motion Field 

The motion field is the vector field induced in the image plane by the relative motion of 
the observer with respect to the environment. The optical flow is the apparent motion of 
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brightness patterns. Under favourable circumstances the optical flow is identical to the 
motion field. The velocity of the image r of a point R is given by 


dr d 1 
— — — “—~R. 
at at R • z 

For convenience, we introduce the notation r t and Rj for the time derivatives of r and 
R, respectively. We then have 


rt 


1 

rTz 


Rt - 


l 

(Rip 


(Rf • z)R, 


which can also be written in the compact form 

r<= (r-'^ ( £>< ( Rt xR ))’ 


since a x (b x c) = (c • a)b — (a • b)c. The vector rt lies in the image plane, and so 
(r t • z) = 0. Further, rt = 0, if Rt || R, as expected. 

Finally, noting that R = (R • z)r, we get 


rt 


j^(zx(R t xr)). 


2.3. Rigid Body Motion 

In the case of the observer moving relative to a rigid environment with translational 
velocity t and rotational velocity u, we find that the motion of a point in the environment 
relative to the observer is given by 

Ri = —w x R — t. 

Since R = (R • z)r, we can write this as 

R( = -(R-z)wxr-t. 

Substituting for Rf in the formula derived above for r*, we obtain 

r « = -(zx (rx (rxw-^rt))). 

It is important to remember that there is an inherent ambiguity here, since the same 
motion field results when distance and the translational velocity are multiplied by an 
arbitrary constant. This can be seen easily from the above equation since the same 
image plane velocity is obtained if one multiplies both R and t by some constant. 


2.4. Brightness Change Equation 

The brightness of the image of a particular patch of a surface depends on many factors. 
It may for example vary with the orientation of the patch. In many cases, however, it 
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remains at least approximately constant as the surface moves in the environment. If we 
assume that the image brightness of a patch remains constant, we have 



dE dr dE_ Q 
dr dt dt ’ 

where dE/dr = ( dE/dx,dE/dy,0) T is the image brightness gradient. It is convenient 
to use the notation E r for this quantity and Et for the time derivative of the brightness. 
Then we can write the brightness change equation in the simple form 


E r • rt + Et = 0. 


Substituting for r t we get 

E t -E e ■ (zx (r x (r xu - = 0. 

Now 

E r • (z x (r x t)) = ( E r x z) • (r x t) = (( E r x z) x r) • t, 
and by similar reasoning 

E r ■ (z x (r x (r x u)) 'j = [((E r xz)xr)xr|-u. 


so we have 

Et — [({Er xz)xr)xr|'W + ((^' r x *) x r ) • t = 0. 

To make this constraint equation more compact, let us define c = Et, s = ( E r x z) x r, 
and v = -8 x r; then, finally, 


c + v • w + ——-s • t = 0. 

R ■ z 

This is the brightness change equation in the case of rigid body motion. 


2.5. Planar Surface 

A particularly impoverished scene is one consisting of a single planar surface. The equa¬ 
tion for such a surface is 

R • n = 1, 

where n/ |n| is a unit normal to the plane, and 1/ |n| is the perpendicular distance of the 
plane from the origin. Since R = (R • z)r, we can write this as 

1 

rn = r, 

R • z 

so the constraint equation becomes 

c + v • u) + (r • n)(s • t) = 0. 
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This is the brightness change equation for a planar surface. 

Note again the inherent ambiguity in the constraint equation. It is satisfied equally 
well by two planes with the same orientation but at different distances provided that the 
translational velocities are in the same proportions. 


3. Recovering Motion and Structure 

Given image brightness E(x,y,t), and its spatial and time derivatives, E e and Et, over 
some region I in the image plane, we are to recover the translational and rotational 
motions, t and w, as well as the plane n. Using the constraint equation developed above, 
we could do this using image information at just a small number of points. At each point 
we get one constraint and we have nine unknowns to recover—or rather, eight, since we 
can recover the distance of the plane and the translational velocity only up to a scale 
factor. 


3.1. Least Squares Approach 


Image brightness values are distorted with sensor noise and quantization error. These 
inaccuracies are further accentuated by methods used for estimating the brightness gra¬ 
dient. Thus it is not advisable to base a method on measurements at just a few points. 
Instead we propose to minimize the error in the brightness constraint equation over the 
whole region I in the image plane. So we wish to minimize 


J = j j [c + v • w + (r • n)(s • t)] 2 dx dy 


by suitable choice of the translational and rotational motion vectors t and u>, as well as 
the normal to the plane, n. We are particularly interested in the question of uniqueness 
of the solution and in computational efficiency. 

For an extremum of J we must have 


dJ dJ 

dw ~ °’ dt~ ’ 


and 


dJ 

dn 


0 . 


That is, 


II 

II 


II 


c + v • w -f (r • n)(s • t)] v dx dy = 0, 


( r ' n ) [c + v • w -f (r • n)(s • t)] adxdy = 0, 
(s • t) [c + v • w + (r • n)(s • t)] r dx dy = 0. 


These equations comprise nine non-linear (scalar) algebraic equation in terms of the 
observer motion, t and w, and the surface normal n. 
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4. Iterative Solution Methods 






To recover the motion vectors t and u, and the surface normal n, we have to solve a 
set of non-linear equations, which we will call the planar motion field equations. Some 
observations about these equations are in order: The first equation is linear in w, t, and 
n. The second equation is linear in u and t, but quadratic in n. Finally, the last equation 
is linear in to and n, but quadratic in t. We will exploit the linearity of these equations 
to formulate two iterative schemes. 

4.1. First Scheme 

We can rearrange the motion and surface recovery equations to get 

^JJ (vx T )dxdy u + ^JJ( r ' n )(vs r ) dx dy t = -JJcvdxdy, 

^JJ (r • n)(sv T ) dxdy u + ^JJ( r ' n) 2 (ss T ) dxdy t = — JJ c(r • n)s dxdy, 

\^JJ (s • t) 2 (rr r ) dx dyj n = — JJ[c + (v • w)] (s • t)r dx dy, 

The first pair can be grouped in the form 

/Mi M 2 
\M% M 4 

where 

Mi = JJ(xx T )dxdy,- M 2 = JJ (r • n)(vs r ) dx dy, M 4 = J J (r • n) 2 (ss T ) dx dy, 
* -II cxdxdy, and d 2 = JJ c(r • n)s dxdy. 

This can be solved for t and w, given the surface normal n. The last equation is 

N 4 n = -g, 

where 


)(:)-- 


dt\ 

d 2r 


g 


N 4 = JJ (s-t) 2 (rr T ) dxdy, 

= JJ [c + (v-w)] (s-t)r dxdy, 


and can be solved for the surface normal n, given the pair of vectors t and w. 
The motion vectors are given by 

u = (Ml 1 M 1 - M 4 _1 M ty 1 (M 4 *d 2 - MJ J di), 

t = (m~ x M 2 - M- r M 4 ) _1 ^Mj r d 2 - M^di) , 
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where ( T ) denotes the inverse of the transpose of a matrix. This can also be written in 
the form 

t = (m 4 Ml - d 2 ) , 

u; = -MJ -1 (di + M 2 t). 

° r ’ _i 

w=(Mi- M 2 M 4 " 1 M l ) (M 2 M ^ M 2 - di), 

t = -MJ 1 (d 2 + M^o;), 

The surface normal is simply given by 

n = -N^g. 

All arrays are either 3x3 matrices or vectors of length 3, and therefore, the solutions for 
u), t, and n can be computed easily. Actually, most of the indicated matrix inversions do 
not have to be carried out explicitly, since it is computationally cheaper to solve these 
linear matrix equations by elimination. 

So in summary: We start with an initial guess for n. Using above equations, we solve 
for t and io in terms of the current value of n, and then for n in terms of the current 
values of t and w. After this, we evaluate the improvement in the solution to either go 
to next iteration or stop if the solution has not improved. 


4.2. Second Scheme 


The first pair of the motion and surface recovery equations depend linearly on t and w. 
As before, 



\JJ (w T ) dx dy 

u + ^JJ \r ■ n)(vs r ) dx dy 

t = — JJ cv dx dy, 

.J 

fj(r ■ n)(sv r ) dx dy 

u + J (r • n) 2 (ss T ) dxdy 

t = — J J c(r • n)s dx dy, 

which can be solved for t and w in terms of n. Furthermore, the first and last equations 
depend linearly on n and w: 


[ff^ T )dxdy 

w + [//(. • t)(vr T ) dxdy 

n = - JJ cv dxdy, 


IJ (s • t)(rv r ) dxdy 

u + [//(. • t) 2 (rr r ) dx dy 

n = - JJ c(s • t)r dx dy. 


Given t, these may be solved for n and u. For simplicity, let Mi, M 2 , M 4 , N 4 , dj, and 
d 2 be as defined earlier, and let: 


Ni = Mi, di=ei, 

N 2 = JJ { s ' t)(vr r ) dx dy, and e 2 = JJ c(s ■t)rdxdy. 
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Then 

(M x M 2 \/ w \ /dA 

U 2 r U/’ 

and 

/Nj Nj\/ W \ (»i 

VNj NjW \ e 2 

The solution of the above equations is given by 

w - (M4M2 -MjMj), 

t = (M^Mi - r M 4 ) 1 (Mj :r d 2 - , 

and 

U = (nj-'N, - N 4 - 1 N 2 1 ')' 1 (Nj‘e 2 - NjV), 
n = (n^Nj - NJ t N 4 )"‘ (NJ r e, - Nj‘ ei ) , 

These may be rewritten in either of two asymmetrical forms as shown above. 

Again, most of the indicated matrix inversions do not have to be carried out explicitly, 
since we can solve the equations by elimination. 

In this scheme: We start with an initial guess for n. We solve for t and u in terms 
of the current value of n, and update t, then solve for n and u in terms of the current 
value of t, and update n, and, finally, evaluate the improvement in the solution to either 
continue with the next iteration or stop if the solution has not improved. 



4.3. Division of Labour 




These methods would not be very attractive, if we had to perform integrations over 
the whole image region I during each iteration, in order to collect the matrices and 
vectors appearing in the equations. Fortunately, this is not neccessary. One can see 
this by writing the equations for the components of the matrices and vectors using the 
summation convention of tensor calculus (that is, there is an implicit summation over 
any index which appears twice in an expression): 


<*>• = // 


{Mi},j = JJviVjdxdy, 

{M 2 },j = \jf v t s J r k dx dy n k , 

= JJi s.Sjrkri dx dyj n k n h 

{d2}<= [//“ 


cvj dxdy, and 


•j dx dy 


n 
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r v 


and 


{Ni} t j = JJjViVj dxdy, 

{N 2 } t j = ^JJ Vis k rjdxdy t k , 

{N 4 },,y = ^JJ^s k sirirjdxdy t k ti, 

{ei}< = JJcvjdxdy, and {e 2 } t = j csyrj dxdyj 

{g}> = III criSj dx dyj tj + | J J s k r{Vj dx dyj t k Wj. 




Mi = Ni and di = ei do not depend on u, t, or n, and so need only be computed once. 
Also, ( cvi ), (tyuy), (cSjTy), ( r k V{Sj ), and (r*.r/ 5 ,sy) depend only on r, E r , and E t) and so 
can be integrated over the image once. This appears to be a set of 3 + 9 + 9 + 27 + 81 = 129 
numbers, but, because of symmetry in (u,i>y), and (r^r/s^sy), only 81 numbers have to be 
stored. These accumulated totals represent all the image information needed to solve the 
motion recovery problem. 

In the first scheme, we only perform 279 multiplications per iteration; The updating 
of the coefficients of the planar motion field equations involves 27 + 9 + 42 + 42 + 42 = 162 
multiplications to compute M 2 , d 2 , M 4 , N 4 , and g (note that M 4 and N 4 are symmetric). 
The updating of w, t, and n, in comparison, requires 117 multiplications. 

In the second scheme, 696 multiplications are carried out at each iteration; We com¬ 
pute the matrices M 2 , M 4 and the vector d 2 , required for the first half of the iteration, 
in 27 + 42 + 9 = 78 multiplications. The same number of multiplications is needed to 
compute the matrices N 2 , N 4 and the vector e 2 required in the second half. Further, 
solving for w and t takes about 270 multiplications, as does solving for u> and n in the 
second half of each iterative step. 

Through a selected example, we will show that the second scheme has a much better 
convergence rate at the expense of more computation per iteration. 


5. Uniqueness 

It is important to establish whether more than one solution is possible. In general, this is 
clearly so, since an image of uniform brightness could correspond to an arbitrary, uniform 
surface moving in an arbitrary way. So the brightness gradients, or lack of brightness 
gradients, can conspire to make the problem highly ambiguous. What we are interested 
in here is whether two different planar surfaces can give rise to the same motion field 
given two different translational and rotational motions of the imaging system. 

In our terms then, the question becomes: given that the brightness change equation 
is satisfied for the motion t and w and the planar surface n, is there another motion t' 
and tv' and another planar surface n' that satisfies the same equation at all points in 
the region I and for all possibly ways of marking the surface? Note that we have to 
consider a whole image region, since the problem is underconstrained if we only have 
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information along a line or at a point in the image. We also have to include the condition 
that the constraint should be satisfied for all possible surface markings to avoid the kind 
of ambiguity discussed above, where brightness gradients fortuitously line up with the 
motion field to create ambiguity. 


5.1. Dual Solution 

Suppose that two motions and two planar surfaces satisfy the brightness change equation. 
Then, we have 

c + v • oj + (r • n)(s • t) = 0, 
c + v • u' + (r • n')(s • t') = 0. 

Subtracting these equations, we get 

v • (w - u/) + (r • n)(s • t) - (r • n')(s • t') = 0. 

Now v — — s x r, so 

—(s x r)• (w — a/) + (r • n)(s • t) — (r • n')(s • t') = 0, 


/“"x 




or 

—r • ((w — a/) x s) + (r ■ n)(s • t) — (r • n')(s • t') = 0. 

If we let u = (wi, then we can write 

( 0 -U >3 +W 2 A 

+W 3 0 —u)\ , 

-U>2 +Wl 0 J 

is a suitable (3x3) skew-symmetric matrix. The element of Cl equals Wk e ijki where 

€ijk, is the permutation symbol. (It equals +1 when the ordered set i, j, k is obtained by 
an even permutation of the set 1, 2, and 3, it equals —1 when the ordered set is obtained 
by an odd permutation, and it is zero if two or more of the indices are equal.) 

Using this notation we can now write 

-r T (n - n')s + r r (nt r )s - r r (n't' T )s = 0, 


or just 


-(f) - Cl') + nt r - n't 


,T 


This is to be true for all points r in the image region I and all possible brightness 
gradients. So 


—(n — fl') + nt r — n't ,T = 0, 


where the zero on the right hand side here represents a 3 x 3 matrix of zeros. Now 
ci T = -n , since fl is skew-symmetric, so taking the transpose of the equation we get 


+(n - D') + tn T - t'n' T = 0. 
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Adding the two equations allows us to eliminate (fl — 17'), and we end up with 

(nt r + tn T ) = (n't' 3 " + t'n' T ). 

The trace (sum of the diagonal elements) of nt T is just (n • t), so we see immediately 

rp 

that (n • t) = (n' • t' ). But the above matrix equation involving the dyadic products of 
n and t as well as of n' and t' is much more constraining. 

Consider the following three possibilities: 

( 1 ) If |n'| = 0 or [t'| = 0 , then their dyadic product is a 3 x 3 matrix of zeros. In this 
case the above equation is satisfied if and only if |n| = 0 or |t| = 0. 

(2) If n' || n, t' || t and ]n'| |t'j = |n| |t|, then the two sums of dyadic products are 
equal and the above equation is satisfied. 

(3) If n' || t, t' || n and |n'| |t'| = |n| |t|, then the two sums of dyadic products are also 
equal and the above equation is satisfied. 

It turns out that there are no other ways to satisfy the equation. This can be shown 
using elementary properties of dyadic products (see Appendix) or by inspection of the six 
components of the above equation (because of symmetry there are only six independent 
components). 

The first case above corresponds to purely rotational motion, because either the trans¬ 
lational motion is zero, or the planar surface is infinitely far away, and the translation 
does not generate a perceptible component of the motion field. The solution is unique 
in this case, because we find (fl — f7') = 0, when we substitute back into the matrix 
equation. (This is nothing new, since it has been known for some time that the solution 
is unique in the case of purely rotational and purely translational motion [3].) 

In the second case we find that nt T = n't' T , since the vectors are parallel and the 
product of their size is constrained by the condition n • t = n' • t', derived earlier. Thus 
once again (f7 — O') = 0. Nothing new is obtained here, since we already know that 
we can change the lengths of the vectors n and t as long as the product of their lengths 
remains constant. 

The third case is the most interesting. Here we have tn r = n't' T so that 

— (f7 - O') + (nt T — tn r ) = 0, 

and thus 

—(fi — n')x + (nt r — tn r )x = 0, 
for an arbitrary vector x. That is, 

x x (w — w') + x x (n x t) = 0 , 
for an arbitrary vector x, so that 

u - w' + n x t = 0, 

or u' = u + n X t. To summarize then, if we ignore scaling of the normal and the 
translational velocity, we obtain a dual solution, given by 

n' = t, t' = n and u 1 = u + n x t. 
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Hay was the first to show the existence of the dual solution [6], although the result has 
apparently been independently rediscovered several times since then [11,13,15]. (The 
most recent of these papers [11], came to our attention only after completion of our 
version of the proof.) 

This dual solution is not different from the original one in the special case that the 
motion is perpendicular to the planar surface, that is, n || t. In this case the solution is 
unique. Further, if t • z = 0, then n' • z = 0. This corresponds to a planar surface parallel 
to the observer’s line of sight, and may be considered to be a degenerate case. 


6. A Selected Example 

We now present the results of a simulation. It is noteworthy to mention that in all 
simulations performed, our algorithms have converged to a solution. However, the number 
of iterations for convergence to a solution depends on the initial condition (as is the case 
with all iterative schemes developed for solving nonlinear equations). In this example, 
we will demonstrate the sensitivity of both schemes to the initial condition. The image 
brightness function was generated using a multiplicative sinusoidal pattern (one that 
varies sinusoidally in both x and y directions), a 45° field of view was assumed, and 
the image brightness gradients were computed analytically to avoid errors due to image 
brightness quantization and finite difference approximations of the brightness gradient. 
In practice, the brightness at image points in two frames would be discretized first, and 
the gradient computed using finite difference methods. 

Table 1 shows the true motion and surface parameters, and the results of a simulation 
that converged to the true solution using the first scheme described earlier. In Table 2, the 
dual solution for the true motion and surface parameters, and the results of a simulation 
that converged to the dual solution are tabulated. In both cases, the solutions after 
various number of iterations are given. The results show that in the first case, the error 
in each parameter after less than 30 iterations is within 10% of the exact value. In 
the second case, this accuracy is achieved in less than 20 iterations. Similar results are 
presented in tables 3 and 4 for the second scheme. Here, very good accuracy is achieved 
in less that 10 iterations for the true solution and about 5 iterations for the dual solution. 

In similar tests, with various motion and surface parameters, accurate results have 
been obtained in less than 40 iterations using the first scheme and a variety of initial 
conditions. The same accuracy for the second scheme required less than 15 iterations. 
More importantly, both schemes eventually converged to one of the two possible solutions. 
However, the results for the particular case where the translational motion vector is 
(almost) parallel to the surface normal have not been as satisfactory. In these cases, 
several hundred iterations were required to achieve reasonable accuracy, even with the 
second scheme. Although the nature of this behavior has not been investigated in detail, 
it appears to resemble that observed when the Newton-Raphson method is applied to a 
problem where two roots are very close to one another. 
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Table 1 : The True Motion and Surface Parameters, and a Summary of the Results of a 
Simulation That Converges to the True Solution. 

True Rotational Motion Parameters w\ = .003 W 2 = .001 103 = —.01 

True Translational Motion Parameters fi = .0005 £2 = —-005 £3 = .0125 

True Parameters of the Surface = .2 ri 2 = .4 713 = 1.0 


Initial Guess for the Simulation ni = 100. ri 2 = 5. ri 3 = — 1. 


Iter. 

(Rotational Par’s) 

(Translational Par’s) 

(Surface Par’s) 


No. 

Wl 

W 2 

W 2 

k 

k 

<3 

ni 

n 2 

713 

10 

.00531 

.00260 

-.01016 

-.00069 

-.00284 

.01301 

.35524 

.1923 

1. 

15 

.00429 

.00178 

-.01008 

-.00006 

-.00384 

.01291 

.27623 

.2742 

1. 

20 

.00353 

.00137 

-.01002 

.00024 

-.00454 

.01270 

.23725 

.3448 

1. 

25 

.00318 

.00117 

-.01000 

.00038 

-.00485 

.01257 

.21718 

.3814 

1. 

30 

.00305 

.00107 

-.01000 

.00045 

-.00495 

.01252 

.20755 

.3945 

1. 

35 

.00302 

.00103 

-.01000 

.00048 

-.00499 

.01250 

.20323 

.3984 

1. 

40 

.00300 

.00101 

-.01000 

.00049 

-.00500 

.01250 

.20137 

.3996 

1. 

45 

.00300 

.00101 

-.01000 

.00050 

-.00500 

.01250 

.20058 

.3999 

1. 

50 

.00300 

.00100 

-.01000 

.00050 

-.00500 

.01250 

.20024 

.4000 

1. 

55 

.00300 

.00100 

-.01000 

.00050 

-.00500 

.01250 

.20010 

.4000 

1. 

60 

.00300 

.00100 

-.01000 

.00050 

-.00500 

.01250 

.20004 

.4000 

1. 

65 

.00300 

.00100 

-.01000 

.00050 

-.00500 

.01250 

.20002 

.4000 

1. 

70 

.00300 

.00100 

-.01000 

.00050 

-.00500 

.01250 

.20000 

.4000 

1. 


Table 2: The Dual Motion and Surface Parameters, and a Summary of the Results of a 
Simulation That Converges to the Dual Solution Using the First Scheme. 

Dual Rotational Motion Parameters wi = .013 W 2 = —.001 W 3 = —.0112 
Dual Translational Motion Parameters ti = .0025 £2 = -005 k = .0125 
Parameters of the Dual Surface ni = .04 ri 2 = —.4 ri 3 = 1.0 


Initial Guess for the Simulation n\ = .5 «2 = 1.5 713 = —1. 


Iter. 

(Rotational Par’s) 

(Translational Par’s) 

(Surface Par’s) 


No. 


W 2 W 3 

k 

k 

k 

ni 

ri2 

”3 

10 

.01302 

-.00120 -.01118 

.00266 

.00503 

.01247 

.01941 

-.4021 

1 . 

15 

.01299 

-.00108 -.01119 

.00256 

.00500 

.01249 

.03220 

-.3992 

1 . 

20 

.01299 

-.00103 -.01120 

.00253 

.00500 

.01250 

.03692 

-.3993 

1 . 

25 

.01300 

-.00101 -.01120 

.00251 

.00500 

.01250 

.03876 

-.3996 

1 . 

30 

.01300 

-.00101 -.01120 

.00250 

.00500 

.01250 

.03950 

-.3998 

1 . 

35 

.01300 

-.00100 -.01120 

.00250 

.00500 

.01250 

.03980 

-.3999 

1 . 

40 

.01300 

-.00100 -.01120 

.00250 

.00500 

.01250 

.03992 

-.4000 

1 . 

45 

.01300 

-.00100 -.01120 

.00250 

.00500 

.01250 

.03997 

-.4000 

1 . 

50 

.01300 

-.00100 -.01120 

.00250 

.00500 

.01250 

.03999 

-.4000 

1 . 
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Table 3: The True Motion and Surface Parameters, and a Summary of the Results of a 
Simulation That Converges to the True Solution Using the Second Scheme. 

True Rotational Motion Parameters = .003 *102 = .001 11)3 = —.01 

True Translational Motion Parameters tj. = .0005 £2 = —-005 £3 = .0125 

True Parameters of the Surface ni = .2 ri 2 = .4 TI 3 = 1.0 


Initial Guess for the Simulation ni = 100. ri 2 = 5. ri 3 = — 1. 


Iter. 

(Rotational Pax’s) 

(Translational Par’s) 

(Surface Par’s) 


No. 


w 2 

Wz 

*1 

h 

h 

ni 

n 2 

^3 

5 

.00254 

.00105 

-.00990 

.00039 

-.00546 

.01202 

.20581 

.44259 

1. 

10 

.00296 

.00100 

-.00999 

.00049 

-.00504 

.01246 

.20039 

.40375 

1. 

15 

.00300 

.00100 

-.01000 

.00050 

-.00500 

.01250 

.20004 

.40037 

1. 

20 

.00300 

.00100 

-.01000 

.00050 

-.00500 

.01250 

.20000 

.40004 

1. 

25 

.00300 

.00100 

-.01000 

.00050 

-.00500 

.01250 

.20000 

.40000 

1. 


Table 4: The Dual Motion and Surface Parameters, and a Summary of the Results of a 
Simulation That Converges to the Dual Solution Using the Second Scheme. 

Dual Rotational Motion Parameters wi = .013 w% = —.001 W 3 = —.0112 
Dual Translational Motion Parameters t\ = .0025 £2 — -005 £3 = .0125 

Parameters of the Dual Surface n\ = .04 ri 2 = — .4 ri 3 = 1,0 


Initial Guess for the Simulation ni = .5 n 2 = 1.5 ri 3 = —1. 


Iter. 

(Rotational Par’s) 

(Translational Par’s) 

(Surface Par’s) 


No. 

W\ 

U>2 

h 

h 

*3 

ni 

n 2 

*3 

5 

.01330 

-.00102 -.01122 

.00248 

.00531 

.01221 

.03790 

-.42663 

1 . 

10 

.01303 

-.00100 -.01120 

.00250 

.00503 

.01248 

.03979 

-.40245 

1 . 

15 

.01300 

-.00100 -.01120 

.00250 

.00500 

.01250 

.03998 

-.40024 

1 . 

20 

.01300 

-.00100 -.01120 

.00250 

.00500 

.01250 

.04000 

-.40002 

1 . 

25 

.01300 

-.00100 -.01120 

.00250 

.00500 

.01250 

.04000 

-.40000 

1 . 



14 


Motion of Planar Objects 


7. Summary 

The problem of recovering the motion of an observer relative to a planar surface directly 
from the changing images (direct passive navigation) was investigated, and formulated 
as one of unconstrained optimization. Using conditions for optimality, it was reduced 
to solving a set of nine simultaneous non-linear equations, which we termed the planar 
motion field equations. Two iterative schemes for solving these equations were presented. 

It was shown that all information in the image concerning motion recovery can be cap¬ 
tured by the moments of the image brightness derivatives which constitute the coefficients 
of the motion and surface recovery equations. These moments are computed during an 
initial pass over the relevant image regions so that there is no need to refer to the image 
after every iteration. This reduces the computation to storing 81 moments and perform¬ 
ing less than 300 multiplications per iteration in the first scheme and approximately 700 
multiplications in the second scheme. 

We also gave a compact proof that the problem can have at most two planar solutions. 
Through a selected example with synthetic data, it was shown that both schemes may 
converge to either of the two solutions, depending on the initial condition. In practice, 
once a solution is obtained, the other can be computed using the equations given for the 
dual solution. 

In the tests carried out, both algorithms have converged to a possible solution, and 
accurate solutions have been obtained in less than 40 iterations using the first scheme, 
and in less than 15 iterations in the second one. As mentioned earlier, the results have 
not been as satisfactory when the translational motion component is perpendicular to 
the planar surface. These cases required several hundred iterations of either scheme for 
accurate solutions. 

Even though both schemes require approximately the same number of computations 
for convergence to a solution (second scheme converges faster but requires more compu¬ 
tation at each iteration, as discussed earlier), the second one seems more appropriate for 
parallel implementation. 
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Appendix 

Lemma 1: (ab T — ba r )c = c x (a x b). 

Proof: (ab r — ba r )c = a(a • b) - b(a • c) = c x (a x b). I 

Lemma 2: If ab T = cd r , then either 

(1) |a| = 0 or |bj = 0, and, |c| — 0 or |d| = 0, or 

(2) a || c and b || d. 

Proof: c r (ab r )d = c r (cd r )d, or (c • a)(b • d) = (c • c)(d • d). Now if |a| = 0 or |b| = 0, 
then |c| 2 |d| 2 = 0, which implies that either |cj = 0 or |dj = 0. Conversely, if |c| = 0 or 
|d| = 0, then |a| = 0 or |b| = 0. 

For the rest of the proof we assume that |a|, |b|, |c|, and |d| are non-zero. Now 
(ab T )b = (cd r )b so that a(b • b) = c(d ■ b). It follows that a || c. Similarly, a 3 ’(ab T ) = 
a r (cd T ) so that (a • a)b r = (a • c)d 3 \ Therefore, b || d. I 


Lemma 3: If M and N are (3 x 3) skew-symmetric matrices and M 2 = N 2 , then 
M = ±N. 


Proof: The result is immediate if M and N are zero. So from now on we assume 
that they are non-zero. If M 2 = N 2 , then M 2 x = N 2 x for all vectors x. Using the 
isomorphism between (3 x 3) skew-symmetrical matrices and vectors we have 

m x (m x x) = n x (n x x), 

where the vectors m and n correspond to the matrices M and N respectively. So 


(m • x)m - (m • m)x = (n • x)n - (n • n)x, 


or 


[mm r — (m • m)/]x = [nn r - (n • n)/]x. 

If this is to be true for all vectors x, we must have, [mm 7, — (m • m)J] = [nn T — (n • n )/]. 
Taking the trace of both sides we get m • m = n • n. Consequently, mm r = nn T . Using 
Lemma 2 , we see that m || n. Taken together with |m| = ]n| this means that m = ±n 
and so M = ±N. B 


Lemma 4: If ab T + ba r = cd T + dc r , then either 

(1) |a| = 0 or |b| = 0, and, |c| = 0 or |d| = 0, or 

(2) a || c and b j| d, or 

(3) a || d and b || c. 

Proof: We first show that either ab r = cd r or ab T = dc r . The proof then follows 
using lemma 2. We also need to show that a • b = c ■ d. This follows from the fact that 
Trace(ab r ) = a • b and Trace(cd :r ) = c • d. Now (ab T + ba 3 ’) 2 = (cd r + dc r ) 2 , or 

(a • b)ab r + (b • b)aa r + (a • ajbb 7 + (a • b)ba r 

= (c ■ d)cd r + (d • d)cc r + (c • cjdd 7 + (c • d)dc r . 
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Further, since ab = cd, we have (a-b)(ab 7 + ba 7 ) = (c • d)(cd 7 + dc 7 ). Subtracting 
twice this amount from the previous equality we get 

-(a • b)ab 7 + (b • b)aa 7 + (a • a)bb 7 - (a • b)ba 7 

= -(c • d)cd 7 + (d • d)cc 7 + (c • c)dd T - (c • d)dc r , 

which reduces to (ab r — ba r ) 2 = (cd r — dc 7 ) 2 . Using the result of lemma 3 now, we get 
(ab r — ba 7 ) = ±(cd 7 — dc 7 ). Adding ab 7 + ba 7 = cd 7 + dc 7 to this equality we find 
that either 2ab 7 = 2cd 7 , or 2ab 7 = 2dc 7 . The conclusion follows using Lemma 2. I 
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